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Abstract 

The inversion problem for rational Bezier curves is addressed by using 
resultant matrices for polynomials expressed in the Bernstein basis. The 
aim of the work is not to construct an inversion formula but finding the 
corresponding value of the parameter for each point on the curve. Since 
sometimes one has only an approximation of that point the use of the singular 
value decomposition, a key tool in numerical linear algebra, is shown to be 
adequate. 
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1. Introduction 

The inversion problem (given a point Pq = {xq, yo) on a curve C parametrized 
by P{t) = {x{t), y(t))), compute the value to of the parameter t corresponding 
to this point Pq) is a classical problem in the area of Computer Aided Ge- 
ometric Design (CAGD) which has a pplications, for instance, in computing 
the intersection points of two curves ( JHoscheck and Lasserl . Il993[ ). 

Two pioneering papers on this subject , which also study the problem of 



implicitization, are ( ISederberg et al.l . Il984j : iGoldman et al.l . ll984j ). the use of 



resultant matrices being a fundamental tool in those papers. A nice recent 



surve y of those and related problems can be found in (ISederberg and Zheng 

2nn2h . 
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Two rece nt alternative approaches to the inversion problem have been 
presented in (JGonzalez-Vega and Rual . l2009l : iBuse and D'Andreal . l2006l ) . 

Our aim in this paper is to solve numerically the inversion probl em for 
the c lass of Bezier curves, a very popular type of curves in CAGD (JFarinl . 
2OO2I ) ■ Since for these curves the parametric equations are expressed in the 
Bernstein basis, it is convenient to construct the resul tant matrices using th e 
Bernstein basis directly, an idea already suggested in (JGoldman et al.l . Il984j ) . 

A resultant matrix for polynomials expressed in the Be rnstein basis is the 



so-ca lled Bernstein-Bezout matrix, which has been used in (JBini and Gemignani 



2OO4J ) for designing a fast algorithm for computing the greatest common 
divisor of two polynomials expressed in the Bernstein basis, and then in 
(JMarco and Martinea . 120071 ) for computing the implicit equation of a Bezier 
curve by means of interpolation. An algorithm for the fast computation 
of the Bezout resultant matrix for polynom ials expressed in the monomial 



basis is introduced in (jChionh et al.l . |2002[ ). and an algorithm for the fast 



computation of the Bernstein-Bezout resultant matrix (which we will use in 
our a pproach to the inversion problem) is presented in (JBini and Gemignanil . 
2004h . 

Another fundamental tool in our approach is the singular value decom- 
position (SVD). The SVD is one of the most valuable tools in numerical 
linear algebra and we will employ it in this work for computing an approxi- 
mation of a basis of the nullspace of a Bernstein-Bezout matrix. The use of 
the SVD is specially appropriate when the matrix is not know exactly, what 
happens in our case when only an approximation of the point Pq = (a^o, Vo) 
is known. It must be observed that this is the usual situation, for example, 
when the inversion problem is utilized to solve problems of curve intersection 
flMarco and Martinej . l2004h . 

Related to this situation (the fact that the po int is usually oii l y clos e to, 
but not exactly on, the curve) in Section 10 of (JGoldman et al.l . Il984j ) the 
problem of the geometric interpretation of inversion f ormulae was bri e fly ad - 
dressed, and this problem was later deeply studied in (IWang and Joel . ll995l ). 

It is important to point out that our approach does not carry out any 
transformation between Bernstein and power (mono mial) bas i s, bec ause this 
conversio n could involve a greater los s of accuracy (JFaroukil . Il99ll ) . In this 
sense, in (JBini and Gemignanil . 120041 ) is indicated that for numerical com- 
putations involving polynomials in Bernstein form it is essential to consider 
algorithms which express all intermediate results using this form only. 

The rest of the paper is organized as follows. In Section 2 some basic 



results on the Bernstein-Bezout matrix and the SVD are presented. In Sec- 
tion 3 our procedure for solving the inversion problem for a Bezier curve 
is introduced, leaving for Section 4 its analysis in the more general case of 
rational curves expressed in the Bernstein basis. Section 5 contains some 
numerical experiments that illustrate the performance of our approach, and 
finally Section 6 is devoted to conclusions. 

2. Some basic results 

In this section we recall some basic results on two fundamental tools that 
we will use in this work: the Bernstein-Bezout matrix and the singular value 
decomposition. 

2.1. The Bernstein-Bezout matrix 

The Bernstein-Bezout matrix is a resultant matrix for polynomials ex- 
pressed in the Bernstein basis. It can be defined in an analogous way to the 
definition of the Bezout resultant matrix fo r polynomials expressed in the 
monomial basis ( Sederberg and Zhena . 120021 ). 



Let Pit) = Po/9i")(t) +pi/3i")(t) + . . . + p^/3t\t) y g(t) = go/3i"^(t) + 
QiPi (t) -!-■■■ + Qnl^n {t) be two polynomials expressed in the Bernstein 
basis 

Bn = {/3j"^(t) = ('') (1 - t)"-¥, ^ = 0, . . . ,n}. 



The Bernstein-Bezout matrix B = (bij) G R"-^"- oi p(t) and q{t) is defined by 



pit)q{s) -pis)q{t) 



t- s 



J2k„pt^'^(t)p^^\s) 



which can be equivalently rewritten as 



p{t)q{s) -pis)q{t) 



t-s 



Pt'\s) ... Pt-,'\s))B 



.(n-l) 



V Pri'it) J 



The following result is a direct consequence of the definition of the Bernstein- 
Bezout matrix: 



Theorem 1. Let p{t) and q{t) be two polynomials expressed in the Bernstein 
basis Bn, and B be the Bernstein- Bezout matrix of p{t) and q(t). If to is a 
common root of p{t) and q{t), then det{B) = 0. 

A fast algorithm for computing the entries of the Bernstein-Bezout ma- 
trix (an algorithm •which w e will use in this work) has been presented in 



( iBini and Gemignanil . |2004| ). Since, for reasons explained in Section 5, we 
will use for the examples of that section the symbolic computation system 
Maple for the computation of the Bernstein-Bezout matrix, we include here 
the the Bini and Gemignani algorithm written in the Maple language: 

for i from 1 to n do 

B[i,l] : = (n/i)*(p[i]*q[0]-p[0]*q[i]); 
od; 

for j from 1 to n-1 do 

B[n,j+1] := (n/(n-j))*(p[n]*q[j]-p[j]*q[n]) 
od; 

for j from 1 to n-1 do 
for i from 1 to n-1 do 

B[i,j+1] : = (n"2/(i*(n-j)))*(p[i]*q[j]-p[j]*q[i]) 
+((j*(n-i))/(i*(n-j)))*B[i+l,j]; 
od; 
od; 

2.2. The Singular Value Decomposition 

Let us now recall the concept of singular value decomposition (SVD): 
Given anmxn real matrix A, there exist orthogonal matrices U (of order 
m) and V (of order n), and a diagonal matrix S (of size m x n) such that 

A = uj:v^. 

This factorization of A is called the singular value decomposition (SVD) of 
A. 

The r (with r < m,n) nonzero diagonal entries of S are the singular 
values of A ( i.e. the positive square roots of the eigenvalues of A^A). If 
there are r (nonzero) singular values, then r is the rank of A. 



But we also obtain a very important additional advantage from the com- 
putation of the SVD of a matrix A of size m x n (including V, not only S) : 
the last n — r column s of of V form a basis of the nullspace of A (see (IStrand . 
19881 : iDemmei Il997f )l 

The SVD provides the best way of estimating the rank of a matrix whose 
entries are floating point numbers, specially in the presence of round-off er- 
rors. This fact was clearly stated in a historical paper by Golub and Kahan 
( jGolub and Kahanl 119651 ). where we can read the following sentences: "In 
the past the conventional way to determine the rank of A was to convert 
A to a row-echelon form... But in floating-point calculations it may not be 
so easy to decide whether some number is effectively zero or not... In other 
words, without looking explicitly at the singular values there seems to be no 
satisfactory way to assign rank to A" . 

Good genera l references in connection w ith th e computation o f the SVD 
of a matrix are ( JGolub and Van Loam Il996l ) and (JDemmell . Il997l ) . 



3. The inversion problem 

Let P{t) = {x{t),y{t)) be a proper parametrization of a plane rational 
Bezier curve C given by: 



x(t) 



yit) 



?(")/ 



An) 



in), 



wpaoPrit) + wiam'jt) + ■■■ + WnanPrit) 

w,P^;\t) + WrPt\t) + ■■■ + WnPt\t) 
wohol3t\t) + wihil3t\t) + ■■■ + Wnbn/3t\t) 



te[o,i]. 

A parametrization P{t) = {x{t),y{t)) of a curve C is said to be proper ii 
every point on C except a finite number of exceptional points is generated 
by exactly one value of the parameter t. Since every rational curve has a 
proper parametrization, we can assume that the considered parametrization 
is proper. Several results on the prop erness of curve parametrizations can be 
found in (ISendra and Winklerl . l200l[ l . 

Our aim in this section is to solve the inversion problem for the case of 
a non exceptional given point Pq = (xo,|/o) G C, where C is a Bezier curve 
properly parametrized by P(t), that is to say, to compute the single value to 
such that -P(to) = Pq- The following result will be essential in our solution 
to this inversion problem. 



Theorem 2. Let p{t) and q{t) he two polynomials expressed in the Bern- 
stein basis Bn- If to is a common root of p{t) and q{t), then the vector 

{(^0 (^o))/^i ~ (to) , ■ ■ ■ , PnSi (^o)) is in the nullspace of the Bernstein- 
Bezout matrix of p{t) and q{t). 

Proof. Let B be the Bernstein-Bezout matrix of p{t) and g(t), B be the 
Bezout matrix of the polynomials p{t) and q{t) expressed in the monomial 
basis {1, t, t^, . . . , t"}, and A^ be the lower triangular matrix of change of 



basis from the Bernstein basis Bn-i to the power basis |1, t, t^ , . 
The following relationship is satisfied (JBini and Gemignanil . |2004| ): 



r-i}. 



B = NBN^. 
Let to be a common root oi p{t) and q{t), then 



/ 1 \ 

to 



B 



,n-l 








(see (jSederberg and Zhengl . |2002| )). In this way 



/ 1 \ 



NBN^ 



and, as A^ is a regular matrix 



/0\ 





V ^0-' / V / 



BN^ 



I ^ \ 




/o\ 


to 







tl 


= 





j.n— 1 
V ^0 / 




^oj 



Therefore the vector 



/ 1 \ 



N^ 



''0 



n-1 



V to ) 

belongs to the nullspace of the Bernstein-Bezout matrix B. Taking into 
account that A^ is the matrix of change of basis from the Bernstein basis 
Bn-\ to the power basis {1, t, t^, . . . , t""^} 



/ 1 \ 

t 



N' 



\ 



t' 



in— 1 






and in consequence the vector (/3q (^o)? hi (to), . . . , /3^_-^ (^o)) is in the 
nullspace of B, the Bernstein-Bezout matrix of p{t) and q{t). 



From now on, we will consider 

p{t) = Woao/3iP\t) + wiai/3^r\t) + ■ 



+ Wnan/3n\t) 



- yoiwoPi;'\t) + w,p[''\t) + ■■■ + wA^\t)), 

and let B be the Bernstein-Bezout matrix of these two polynomials p{t) and 
lit)- 

Taking into account Theorem 2, the vector {/3q (to),/?}" (^o); 
1^2 (^o), ■ ■ ■ , Pn-i (to)) is in the nullspace of B. Moreover, as Pq = (xq, yo) 
is a non exceptional point of C, which is properly parametrized by P{t), the 
rank of -B is n— 1 and therefore the vector {I3q (to), /3|"' (to); • • • ? Pn-i (to)) 
is a basis of the nullspace of B. 

Let us observe that it is very usual not to know the point Pq exactly. In 
this way, considering the results included in Section 2.2, it is convenient to 



7 



compute an approximation of a basis of the nullspace of B by using the SVD. 
Naturally, we do not have the exact matrix B but an approximation of B 
that we will denote hy B. A theoretical study of the application of SVD to 
the computation of an ap proximation of the nullspace can be seen in Section 
1 of Chapter 5 of lStrand {1988). 

As the matrix V of the SVD of B is orthogonal, the last column of V gives 
us an approximation of a multiple of the vector (/Sg ~ (to),/3i ~ (to), ■ ■ ■ ^ l^n-i (^o)) 
with Euclidean norm equal to 1: 



Since 



[zo,zi, . . . ,Zn-i) — a{f3Q (to), (3-^ (to), • • • ,/3n-i (^o))- 



z^ Pt'\t) ("7^)(l-t)"-^-^f {n-^)t 



we obtain the value of the parameter t we are interested in computing: 

izi 



t, 







iZi + (n- i)Zi-i 



1,2,. ..,n 



1. 



4. A more general situation 



In this section we show how our procedure can also be used for the case of 
a non Bezier rational curve C properly parametrized by P(t) = {x{t),y{t)), 
where x{t) = ^, y{t) = ^ and Ui{t) {i = 1,...,4) are polynomials 
expressed in the Bernstein basis o„. In this situation some care must be taken 
when deg{p) ^ deg{q), where p(t) = ui(t) — xu2{t) and q{t) = us(t) —yuiit). 

Given a polynomial pit) expressed in the Bernstein basis Bn, we denote 
by deg{p) the degree of p{t) expres sed in the power basis. It must be noticed 
that n > degip) (see Section 3 in ( JBuse and Goldmanl . l2008l )). 



The following example (JMarco and Martinezl . 120071 ) serves to show what 
happens in this special situation: 

Example 1. Let B4, be the Bernstein basis of the space of polynomials 
of degree less than or equal to 4, and let us consider the curve given by the 
parametric equations 



x{t) 



?(4)^ 



?{4) 



^(4)/ 



(4)/ 



Ap'^\t) + Apr it) + 3/3r(t) + m\t) + writ) 



?(4)^ 



?(4)/ 



?(4)/ 



?{4) 



^rit) + pt\t) + /3r(t) + /3r(t) + 3/3r(t) 



(4)/ 



?(4)/ 
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In this case 

p(t) = (4-x)/3j')(t) + (4-x)/3;')(t) + (3-x)/3j')(t) + (3-x)/3f(t) + (7-3x)tf)(t) 
and 

However, if we write p and q in the power basis we have 

p{t) =A-x-6t^ + 8t^ + {-2x + l)t^ 

(a polynomial of degree 4 in t), while 

q{t) = 2-y + At-6f + 4t^ 

a polynomial of degree 3 in t. In this example, our approach works per- 
fectly when computing the parameter to corresponding to every point Pq = 
{xo,yo) G C such that xq ^ ^, i.e., it works for every point Pq E C except 
for the point Pq = (1/2,-3.0395517). The reason why we have problems 
when xo = I is that xq = | makes the leading coefficient of p{t), and 
the nuUspace of the corresponding Bernstein-Bezout matrix has dimension 
greater than 1. When this happens, a possible solution is to consider the 
Sylvester matrix of the polynomials p{t) (expressed in the Bernstein basis 



Bj) a nd q{t) (expressed in the Bernstein basis B3) ( IWinkler and Goldmanl . 



20031 ) instead of the Bernstein-Bezout matrix oip{t) and q{t) expressed in the 
Bernstein basis B4, and then proceed in the same way as we have described 
for the Bernstein-Bezout matrix. 

5. Numerical examples 

In this section we include several examples illustrating the performance 
of our procedure. A detailed description of the way in which we proceed at 
each step of our approach is included in the first one. 

Example 2. Let us consider a rational Bezier curve C given by 



E;:„».(i)**(i-*)"- 



where the points of the control polygon are 

(ao,6o) = (14,14), (ai,6i) = (11,15), (02,62) = (9, 15), (03,63) = (7, 15), 

(04,64) = (4,14), (05,65) = (3,12), (O6,66) = (3,10), (07,67) = (7,8), 

(og, bs) = (4, 6), (09, 69) = (14, 4), (010, 610) = (12, 2), (on, 611) = (8, 2), 

(ai2, 612) = (6, 2), (013, 613) = (4, 3), (014, 614) = (3, 4), (015, 615) = (2, 5), 

and the weights are 

Wo = 2, wi = 2, W2 = 2, W3 = 1, W4 = 2, W5 = 5, wq = 5, wj = 1, 

W8 = 3, Wg = 3, Wio = 3, Wu = 3, W12 = 2, W13 = 1, Wu = 1, Wi5 = 1. 

Let us consider Pq = (8.50665, 14.3420), a point in C that we know approxi- 
mately (it is an approximation of the point corresponding to the value of the 
parameter t = y). We take Pq = {xq = ^^^^,yo = ^^^) and we construct 
the Bernstein-Bezout matrix oip{t) and q{t) by means of the Maple code in- 
troduced in Section 2.1 using exact arithmetic. The reason way we use exact 
arithmetic instead of fl oating point arithmet i c is th e following: although the 



algorithm presented in ( iBini and Gemignanil . |2004| ) for computing the entries 
of the Bernstein-Bezout matrix is fast (it has a computational cost of 0(r2^) 
arithmetic operations) this algorithm does not have high relative accuracy, 
and therefore it is possible that important errors appear in the computation 
of some entries of the matrix. 

Once we have in Maple the approximation B of the Bernstein-Bezout 
matrix oi p{t) and q{t), we put it in Matlab where we compute its singular 
value decomposition. We include here the last column of the matrix V, which 
give us an approximation {(^o; ^1, • • • , ^n-i)} of the basis of the nullspace of 
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B: 

( -2.473682899590338e - 001 \ 

-5.771889609913881e - 001 
-6.253061379846553e - 001 
-4.169368183026464e - 001 
-1.910461999315792e - 001 
-6.368830764276734e - 002 
-1.592182997982856e - 002 
V{:, 15) = -3.033481337354997e - 003 
-4.422976092589377e - 004 
-4.912824640160427e - 005 
-4.044767281991494e - 006 
-2.444466877812876e - 007 
-1.184192736731686e - 008 
-5.593200217925554e - Oil 
\ -4.503186536883329e-011 / 

In the computation of to we use two consecutive components of the vector 
above chosen in the following way: we select the vector component which has 
the greatest absolute value, and then we take among the component before 
and after this the one that have the greatest absolute value. We proceed in 
this way because the components with greatest a bsolute value are usu ally 
the ones which are less sensitive to perturbations (JDemmel et al.l . 120081 ). In 
this case we obtain 

2^2 



tn 



1.428606867264249e - 001. 



1Z2 + 13;zi 

When we compare the value of to we have obtained with the exact value of 
the parameter t = ^ of the exact point 

'78193109744768 131831466405881^ 



A^).y{^) 



eC, 



9191995131007 9191995131007 
corresponding to the given point Pq = (8.50665, 14.3420), we get that the 
relative error we made when using our procedure is 

|l/7-to| 



1/7 



2.48e - 005. 



Taking into account that we have started from an approximation Pq of the 
point P E C with only five exact digits, we can assert that the relative error 
that we have obtained by using our approach is small. 
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Now we present a detailed description of the computations for the non 



Bezier curve presented in Example 1. 



Example 3. Let us consider the curve C introduced in Example 1 and a 
point Pq = (xo, Vo) e C such that xq ^ |, for instance Pq = (3.5542169, 2.8148148) 
which is an approximation of the point corresponding to the value of the pa- 
rameter to = 1/3. We show how our procedure works for these kind of points 
(all the points of C except one). 

If we consider Pq = (^^^^^? '^^^^^)-, the singular values corresponding 
to the matrix B are: 

(71 = 4.212191730287018e + 000, 
02 = 2.075444341475023e + 000, 
as = 5.981428444978487e - 001, 
^4 = 3.357757839963324e - 008. 

As only one singular value (0-4 = 3.357757839963324e — 008) can be consid- 
ered 0, the nullspace of B has dimension 1 and therefore the last column of 
the matrix V of the SVD of B 

( 5.111012703997072e-001 \ 

7.666518828223443e - 001 

3.833259078513024e - 001 
\ 6.388763832491218e - 002 / 



\/(:,4) 



gives us an approximation of the nullspace of the exact Bernstein-Bezout 
matrix B. Proceeding in the same way as in the previous example we obtain 
that 

to = — = 3.333333267311144e - 001, 

a good approximation of t = | if we take into account that Pq is an approx- 
imation of the point corresponding to the value of the parameter to = 1/3 
with only 7 exact digits. 

Now we show what happens when we consider Pq = (1/2,-3.0395517). 
In this case the singular values corresponding to the approximate Bernstein- 
Bezout matrix B are 

o^ = 7.650996902942649e + 001, 
(72 = 4.412324130498519e + 001, 
ag = 8.892019775556218e - 009, 
a4 = 1.013508461498068e - 015. 

12 



Since the singular values a^ and (74 can be considered 0, the nullspace of 
B has dimension 2, and therefore it is not possible to proceed by using the 
approach we have described in Section 3. As we have pointed out in Section 
4, a possible solution for this situation could be to consider the Sylvester 
matrix of the polynomials p{t) (expres sed in the Bernstein bas i s B4.) and q{t) 



(expressed in the Bernstein basis B2) (jWinkler and Goldmanl . |2003[ ) instead 



of the Bernstein-Bezout matrix of p{t) and q{t) expressed in the Bernstein 
basis S4, and then proceed in the same way as we have described for the 
Bernstein-Bezout matrix. 

Finally we present an example with a Bezier curve of small degree. 

Example 4. Let us consider the rational Bezier curve C of degree 3 
given by 

where the points of the control polygon are 

(ao,M = (l,9), (ai,6i) = (2,l), (02, ^2) = (5, 1), (ag, 63) = (4, 1), 
and the weights are 

Wo = 1, Wi = 2, W2 = 2, W3 = 1. 

In this case, as the considered curve has small degree, it is easy to obtain (by 
using the Maple command f f gausselim) the following explicit expression for 
an inversion formula: 

81408X + 13824|/x - 7680x2 _ 276552 - 1008?/ + 1080^^ 



fix,y) 



-216432 - 21024?/ + 67200x + 12672|/x - 6144^^ - 2160?/2' 



The exact point in C corresponding to the parameter t = |isP=(|,^). 
When we consider its approximation Pq = (2.66667, 2.42222) the value of to 
that we get by means of our approach is 

3.333339104290224e-001, 

while the value of to we obtain by using in Maple the above inversion formula 
f{x,y) for Pq is 

3.333319852915495e - 001. 

These results show that, even in the cases we have an explicit inversion 
formula, the computation of the parameter to by evaluating this function 
does not necessarily give better results than using the SVD. 
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6. Conclusions 

In this work we liave sliown liow tlie combination of classical algebraic 
geometry tools (resultant matrices) and numerical linear algebra tools (the 
singular value decomposition) is adequate for addressing the inversion prob- 
lem for a rational Bezier curve, although the more general situation of ratio- 
nal non Bezier curves defined by polynomials in the Bernstein basis is also 
considered. The resultant matrices being used avoid the conversion between 
the Bernstein basis and the monomial basis. 

It must be observed that the aim of our work was not the computation 
of an inversion formula, which is not an easy task for curves of high degree 
and must be approached by using different techniques. Instead we have con- 
sidered the numerical problem of finding the parameter value corresponding 
to a given point on the curve, a point that usually in practice is not known 
exactly. 
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